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I. 



INTRODUCTION 



The parabolic partial differential equation 

-|^ = V* (a(x, t,u) Vu) + f(x,t,u,Vu) (x,t) e J2x(0,t] 

has been the subject of the development of several niimerical 
approximation techniques. The equation is known as the 
heat equation in a specialized form and is useful in the 
study of heat transfer and thermodynamics in general. It 
also has applications in other fields of engineering and 
science, such as fluid dynamics and meteorology. 

The development of Galerkin methods for the solution 
of a general parabolic and hyperbolic equation has been 
relatively recent [3,6,9]. This thesis describes a program 
for the solution of the general parabolic problem using the 
central different Laplace modified alternating direction 
Galerkin scheme described in [5,6,8]. The program was built 
in steps beginning with the one dimensional case and then 
extended to the two dimensional case. 

The program solves the problem for a rectangular domain 
and uses nine Gauss quadrature points for integration in 
the plane. It also allows for nonlinearities in the forcing 
function f(x,t,u,Vu) and the function a(x,t,u) can be 
extended to have nonlinear terms of Vu. The theoretical 
proofs of stability and convergence are available for these 
cases in [4] . 
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The program was written in the Fortran IV computer 
language and tested on an IBM 360/67 computer system. 
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II. 



DEFINITIONS, NOTATION, AND PRELIMINARY CONCEPTS 



In the following section several concepts and definitions 
which are essential for the understanding of the thesis are 
presented. Two specific spaces of functions will be consid- 
ered and definitions of the norms in the spaces will be 
given. 

Let with smooth boundary be the domain of 

interest for the function spaces. The first of these spaces 

2 

is that of square integrable or measurable functions, L (fi). 

2 

If w eL (J2) , then the norm of w is defined as 






/ W^ dx < oo 



2 

L (fi) is a Hilbert space with respect to the inner product 



<f,g> = / f(x) g(x) dx 



0 2 
The Sobolev space H (fl) is equivalent to L (fl) . 



The second space considered will be the Sobolev space 
defined as 



= {w|dPw e (fi) , |p| = 0,1} 



where 




!p 1 = (p^) 



n 





i=l 
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The norm of w is defined as 



n 



II 8x^ II l 2 (fi) ^ 



1/2 



Of special interest for the development of this paper 
1 1 

is the space a subspace of H defined as 



= {w|w eH^(fi) and w = 0, for all x e 9fi} 



The norm for this subspace is defined to be 






n 

Z 

i=l 



3w 

ax 



l 2 (fi) 



1/2 



A major concept which is used in the formulation of 
Galerkin's method is that of the weak solution to a differ- 
ential equation. Let L be the differential operator 
describing the differential equation, that is 



Lu = f 



( 2 . 1 ) 



A function u is said to satisfy the weak form of equation 
(2.1), if for every test function v in some test space 

<Lu,v> = <f,v> (2.2) 

If u is not sufficiently differentiable to be substituted 
in (2.1), but does satisfy (2.2), then u is said to be a 
weak solution. 
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To illustrate the weak solution, consider 



Lu = - u^^ = f(x,t) , xel = (0,1), 

u(0,t) = u(l,t) = 0 ; u(x,0) = Uq(x) , 



t > 0. 

xel 

(2.3) 



The weak form of this equation is 

<u. - u ,v> = <f,v> for all veHj!:(I) (2.4) 

u(0,t) = u(l,t) = 0 
u(x, 0) = Uq (x) 



Integrating by parts results in 






<u , V 
x' X 



> = <f,v> 



for all V 



e hJ(I) 



u(x,0) = Uq(x) 

Note that u £11^(1) can be a solution to 

is a weak solution of (2.3), however a strong 

2 

(2.3) must satisfy ueC (I) for each t. 



(2.5) 



(2.5) , that 
solution of 
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III. 



THE FORMULATION OF THE GALERKIN METHOD 



The formulation of Galerkin methods for the solution of 
differential equations is based on the weak form of the 
equation. Consider the nonlinear parabolic equation 



8u 

3t 



= V* (a(x, t,u) Vu) + f(x,t,u,Vu) , (x,t) efix(0,T] 



u (x, 0) = Uq (x) , X e 



(3.1) 



u(x,t) = 0 (x,t) e3f2x(0,T] 



and the test space Clearly for any veH^Cfi) 



3u 

3t 



V = V • (a (x, t , u) Vu) V + f(x,t,u,Vu)v 



Taking the integral over the region Q yields 



/ V dx = / V* (a(x,t,u) Vu) V dx + / f(x,t,u,Vu)v dx 



Using Green's first formula for integration by parts in 
multiple dimensions results in 



/ V dx = / a{x,t,u)Vu V dx - / a (x, t,u) Vu* Vv dx 

a a 



+ f f(x,t,u,Vu)v dx 
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Since v = 0 on it follows that 



/ a(x,t,u)Vu V dx - 0 

and 

/ V dx + / a(x, t,u) Vu’ Vv dx = / f(x,t,u,Vu)v dx 

ft 

Therefore u satisfies 

,v> + <a (x, t , u) Vu, Vv> = <f (x, t,u, Vu) ,v> (3.2) 

<u(x,0) ,v> = <Uq,v> for all veH^Cft) . 

Let M be a finite dimensional subspace of njcft) . The 
continuous-time Galerkin method is to find for each t e (0,T] 
a differentiable map U ( • , t) : [0 ,T] ^ M , such that 

<1^ ,V> + <a(x,t,U) vu, VV> = <f (x,t,U, VU) ,V> (3.3) 

<U,V> = <Uq,V> , t = 0 , V e M 

Now let M be a subspace of HQ(ft) such that 
M = Span (w^ ,W 2 / . . . /W^) , where ^ linearly 

independent set. Under these conditions it is possible to 
write 

n 

U (x, t) = L a. (t) w . (x) 
j=l J ^ 
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Svibstituting in equation (3.3) and letting V = 

n n 

< L a . ( t) w . (x) ,w. (x) > + <a(x,t,U) I a . (t) Vw . (x) , Vw . (x) > 
j=l J J ^ j=l J J 



= <f (x,t,U,VU) ,w^(x) > 

and 



n 

< I a. (0)w. (x)> = <u„,w. (x)> 
j=l 3 3 u 1 



Since integration is a linear operation and the inner product 
defined is not over the time space, it follows that, 

n ^ n 

E a!^ (t) <w . (x) ,w. (x) > + E a . (t) <a (x, t ,U) Vw . (x) , Vw • (x) > 

j=l J J ^ j=l ^ J ^ 

=<f(x,t,U,VU),w^(x)> (3.4a) 



and 



n 

E a . ( 0) <w . (x) ,w. (x) > = <u^,w.(x)> i=l,2,...,n (3.4b) 

j=l :) J 1 u 1 



Equations (3.4a) and (3.4b) can be written as a system 
of n-nonlinear ordinary differential equations 



Ma' (t) + S(a)a(t) = F 
Ma(0) = b 



(3.5) 
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where S(a) 





= <Wj (x) ,w^ (x) > 



and 



F = (f^) ; = <f (x,t,U,VU) ,w^(x) > 



a = 




T 



b = (<Uq,Wj^>,<Uq,W2> 




T 



Using the properties of the functions w^, it can be 
shown that the matrix M is positive definite, and with the 
additional property that a(x,t,u) is bounded above and below, 
the matrix S can also be shown to be positive definite. If 
in addition a(x,t,u) and f(x,t,u,Vu) satisfy Lipschitz 
conditions, it can be shown that the system of differential 
equations has a \anique solution. 

In the one dimensional case, an approximate solution to 
equation (3.5a) can be obtained by combining it with a 
finite difference scheme to discretize in time as follows. 
Assuming, for simplicity, that a(x,t,u) = a = constant 
and f(x,t,u,u^) = f(x,t) 




m+1 , m 
; + a 



2 



F 



,m+% 



(3.6) 
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This is a system of n equations and n unknowns and can be 
simplified to 



(2M+AtS) = (2M-AtS)a"^ + 2AtF"''‘‘^ 



or 

(2M+AtS) - a^) = -2AtSa’^ + 2AtF"’‘‘’^ 

where o'? = a. (t ) and = <f(t , 7 )/W.> 

In the above equation the time variable t has been discretized 
to tj^ = mAt , where t = T/N , N a positive integer. 

This method is known as the Crank-Nicholson-Galerkin method. 

The one dimensional program was written to set up the 
matrices defined above and solve the algebraic equation (3.7). 
A minor alteration in the program v;ould allow the use of 
(3.7a), and give the advantage of reducing the computational 
roundoff error. The program is also limited to the linearity 
simplifications made above. 

The required accuracy of the method in the two dimensional 
case leads to a large number of basis functions for the space 
and consequently would require excessively large matrices. 

The matrices have a band structure and the band width of the 
matrices is dependent on the number of intervals used in the 
grid set up for the region. These large matrices can be 
avoided by implementing an alternating direction method, and 
using a tensor product basis for the space. This tensor 
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product basis reduces the problem to working with matrices 

which have a constant band width of seven, in the case of 

M being Hermite cubic splines. It also greatly simplifies 

2 

the problem when the region fi R is rectangular. 

Let M = Span(^, (x) , (x) , . . . (x) ) be a basis for 
X 1 z n 

X 

one dimension of the space M and M = Span(n-i (y) /H- (y) , • . . 

y 1 ^ 

ri (y) ) be a basis for the other dimensions of M. M and 
M are bases for one dimensional subspaces of Hj:(I ) and 

j u X 

such that M = ® = Span ( / • • • / , 

...,f X] ) . Also define 

n n 
X y 



<f,g> = / fg dx 

^ I 

X 



<f,g> = / fg dy 



Since M SM forms a basis for the space M, we can write 
X y t- f 



^x ^y 

U(x,y,t) = E Z a (t) (C^(x) Q ri„(y) ) . 



P=1 q=l 



pq 



There are several time discrete m.ethods with which to 
proceed at this point. The program has been designed to use 
the alternating direction version of the centered difference 
Laplace modified Galerkin method described in [5,6,8] and 
defined as 
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,m+l ,,ni-l 

< ,V> + <a"'(x,t,U) VU"', W> 



+ 2 X At < 



2 2 

3 ^2,,m 9 V . ,, ,, 

9 U , = <f (x,t,U,VU),V> 



3x 3y 



9x ay 



ctnd 



<U,V> = <Uq,V> V e My 



(3.7) 



where the superscript of the functions denotes the time 
step at which it is to be evaluated, and 






2u"' + ^ 



and 






X is subject to the restriction 



X > -T- a = -r max a(x,t,u(x,t)) 

X 

t>0 

for stability [5,6]. 

Equation (3.7) can be written in the algebraic tensor 
product form 

(C +2XAtA ) a (C +2XAtA ) = 2(2XAt(C a A +A ac ) 

X X y y X y y y 

+ 4X^(At)^(A^aAy) )a"' 

+ (C^aCy - 2XAt(Cj^aAy + A^aCy) - 4x2(At)^(A^aAy))a"^"^ 
+ 2At<J)"' 
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where 



<I>pq = -<a"'(x,t,U)Vu"',V(Cp nq)> + t,U, VU) , ?p riq> 



and 



C = <C . , C . > 

X X j X 



c, = <n • / h .> 
y 1 D y 



A = 

X X 



A = <riwn^> 

y 1 D y 



Define e'^ = ^ , then equation (3.7) can be 



written as 



(C +2XAtA^) Q (C +2XAtA„) (e"^‘^^-e’^) = -2 (C^ Q C ) + 2At$^ 

X X. y y ^ j 



(3.8a) 



and 



m+1 _ / m+1 . m 

a = (e -e ) + e + a 



(3.8b) 



The scheme (3.8) adds significant efficiency to the program 

and also reduces roundoff error. 

The Laplace centered difference method initially requires 

two sets of the coefficients, those at time zero and at one 

time step later. At time zero the set of coefficients are 

2 

obtained by an L (fi) projection of the initial conditions 
into the space, but the second set presents a problem where 
a variety of methods are available. Two of these methods 
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I 



were attempted, expanding the initial conditions in a 
Taylor's Series about t = 0 and using the Laplace forv/ard 
difference method. The former involved taking the deriva- 
tives of the initial condition function and presented a 
problem when these derivatives resulted in zero. The method 
which was finally used was the Laplace forward difference 
scheme. It is defined as follows 



At first glance the matrices in (3.8a) and (3.9) seem to be 
different, however in the case of this method X is under 
the restriction 



thus a prudent choice of X results in identical matrices and 
it is not necessary to form a second set of coefficient 



the same form as in the previous method and the subroutine 



to the final choice. 

The tensor product of matrices has the property that 



(C^ + XAtA^) S (Cy + XAtAy) - a^) = (At) 



(3.9) 




= max a (x, t , u) 



X e , t > 0 



matrices in the program. It is also noteworthy that has 



used to compute may be used for both with a correction of 
i for the forv/ard difference scheme. These advantages led 



(B QD) = (B fil) (I QD) 
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This property enables a significant simplification of two 
aspects of the programming problem. First, the tensor 
product multiplication on the right hand side of the equation 
can be accomplished by computing first 

(I S Cy) e = Y 

and then 

(C^Q I) y"' 

By properly ordering the vector e'^, this results in a series 

of multiplication by the relatively small matrices and C^, 

that is, matrices which have a constant band width of seven 

for Hermite cubics rather than a band width which depends 

on the number of intervals in the grid of the problem as in 

CSC. 

X y 

This same idea 'can be used in solving the linear system. 
Defining the right hand side of the equations (3. 8, 3. 9) 
to be 8^, let 

(I Q (Cy + 2AAtA^) ) - e^) = Y^ 

Using this, solve first 

((C +2XAtA^) QI)y^ = B"' 

X X 

and then solve 

(I Q (Cy + 2AAtAy) ) - e"') = Y^ 
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Again proper ordering of the vectors involved results in 
the simplification of solving a system of one dimensional 
problems rather than the larger problem of solving the full 
tensor product matrix. 

The two dimensional program has been designed to solve 
the problem in the manner described above. The basis 
functions used in the program are Hermite bicubic splines, 
taken as the tensor products of one dimensional Hermite 
cubic splines, and will be described in Section IV. The 
procedure described above has the drawback that it is in 
general good only on rectangular polygons, and must be 
altered significantly to handle regions other than simple 
rectangles. It can be changed to handle regions which are 
unions of rectangles [5] . 



IV. HERMITE BICUBIC SPLINES 



The Galerkin procedure is dependent on a set of linearly 
independent basis functions for the space considered in the 
problem. Hermite bicubic splines provide such a basis. 

These splines are a high powered tool for the approximation 
of smooth functions. 

The Hermite bicubic splines are taken as the tensor 
product of one dimensional Hermite cubics. These cubics 
have a piecewise polynomial nature, where the polynomial is 
of degree less than or equal to three. Hermite cubics have 
globally continuous first derivatives. Thus, the tensor 
product of the one dimensional cubics has the property that 
both first partials and the first mixed partial derivatives 
are continuous. These properties are sufficient for the 
approximations desired as a result of the program [1] . 

The package of subroutines described in [2] was used for 
the program. It consists of subroutines which compute the 
values of the B-splines at specific points, and given a set 
of coefficients, computes the value of the approximated 
function at specified points. The high degree of flexibility 
of the package was the major reason for its use in the 
program. The program can be extended to use higher degree 
B-splines which gives a greater degree of smoothness and 
accuracy to the approximating function. Also, the extension 
of the program to higher dimensions can be accomplished by 
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taking a higher degree tensor product, and making the 

appropriate alterations in the computational procedures. 

Consider the interval [0,1], and the partition of the 

interval A: 0 = x,v < x, < . . . < x„ = 1 . Let h- = x- -x. ,, 

0—1— — n X xx-1 

and I. = [x. ,,x.]. Define h = max h. , and h = min h. . 

X X- l x j^x' — ^x 

The Hermite cubics are piecewise cubics over for all i, 
and are members of C^(I). The order of accuracy for Ilermite 
cubics is given by the following well-known theorem; 

Theorem . If u eH^(I) and h/h ^ a < °° , then 

inf 11^ " ^11 hS(i) 1 Ch^~^ II ^11 (I) ' 0 1 s f 3 

where M is the space of Hermite cubics. 
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V. ERROR ESTiriATES 



The usefulness of a procedure to obtain an approximate 
solution to a differential equation is measured by its 
degree of accuracy and the efficiency and costs of its 
implementation. The purpose of this section is to derive 
an error bound for the continuous time Galerkin method for 
a linear version of equation (3.1). The following theorem 
and its proof is found in [8]. It is restricted to the case 
where a(x,t) is not a function of u. Proofs for the nonlinear 
case are also available [3,6,9]. 

Theorem 1. Let u be the solution to 



,v> + <a (x, t) Vu, Vv> = 0 veH^Cft) , te (0,T] 



(5.1) 



<u(x,0)>= <Uq(x)> xefi 



and let U be the solution to 



/V> + <a(x,t) VU,VV> = 0 
d t 



V e M , t e (0,T] (5.2) 



<U,V> = <Uq,V> 



t = 0 



and 




(x,t) £ c^ 



(5.3) 
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Then there exists constants C and d, which depend on 
Cq/ c-^, and T, such that 



max 

0<t<T 



I|u-u|Il 2(„) (t) 



+ 



9 J' llu-ull % (t) dt 
0 “o 



< C( 



max 

0<t<T 



u~u 



l 2 (fi) 



T 

+ / ||u-u 

0 






dt 



T 

+ / 

0 



It (“-S>ll^2,n,(t) dti 



(5.4) 



where u is an arbitrary map of (0,T] into W. 

Proof; Replace v with V in (5.1) and subtract (5.2) from 
(5.1) 

< ,V> + <a(x,t) Vu, VV> = 0 

-< 1^ ,V> + <a(x,t)VU,VV> = 0 



results in 



< ,V> + <a(x,t) V(u-U) , VV> = 0 . 

O t 

Let e = u-U and let V = e+e , where e = u-u , then 



< ,e+e> + <a (x, t) Ve, Ve+e> = 0 

o 
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I 



Note 



8e , 



8e 

at 


e dx 


d 

dt 


/ e2 


d 


llell^ 


dt 



l 2 (fi) 



and 



<aVe,Ve> > ||e|| 



from (5.3) 



and 



l<aVe,Ve>l < HaVeH ^ 2 ^^) • ||Veli 



- ^1 11® II hJcq) 



®11 Hj(fi) 



from Schwarz's inequality. 

Young's inequality states |lfH • ||gH < e |lf|| ^ 11<?11^- 

Using this inequality we get 



<ave,vs>| < ||e|| ^ ||s|| 2 



Thus 



Id 1 1 ^ I I ^ a. ^ 3 ^ X X 4- r* 

2 ^ ll®ll l 2{Q) ■*■ *^9t ^ ®( 



®llHj(fi) 



1 ^®1 ll®llnj(fi) ^ 4e I1®I1hJ(J2) 



t e (0,T] 
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and 



2 dt ll®ll ^3t 

- 47 ll®ll nj(fi) ' ^ ^ 

2 

Let c = Cq/2c 2^ , then for C = /c^ 

ll^ll hw * 2 <|| ,6> + C(, l|e|lHl(n) 

- ^ ll®ll hJ(!)) 



Integrating from 0 to t 



L2(fi) 



(t) - II ® II l2 (Qj (0) + Cg / ll®llHQ(f^) 



and 



“ ^ 0 nj(fi) ^ Q ^at 



/ <-||- ,e> dt = <e,e> 

0 



f <e, -||> dt 
0 



|<e,e>| _< el|e|| t 2 /o) ^ ll®ll t .2 



(fi) 
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and 



0 



9 s 1 
9 tl 



dt < 



J 

0 



L 2 (fi) 



(t) dt + 



9s I 
9t 






(t) dt) 



It follows that 



II® II l 2 (fi) " ll®ll L2(ft) ^0 ^II®IIhJ(.Q) 

- Il®ll L2(fi) ■*■ ll®II L^(Q) ^ II®II Hl(f2) 

Q II®IIl2(J^) II® II l 2 (fi) ■*■ II® II l 2 (fi) 

^ f ll|fll^(^2) 



Thus for e sufficiently small 



II 



1 l2(s) *■'* ^ * 0 nj(fi) 

= ll®ll^(a) * ll^ll^(a)<“' ^ f ll«llHi(f!) 
ll^ll L^lf!) <^> ■" ll®ll L2(f!) ^ lllfll l2(!J) 



Gronwall*s inequality : Let f,g, and h be piecewise continuous 

nonnegative functions defined on an interval a < t < b and 
assume that for each t e [a,b] 



f(t) + h(t) < g(t) + / f(s) ds 

0 
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Then 



f(t) + h(t) £ / ^g(s) ds + g(t) 

a 

and if g is nondecreasing 

f(t) + h(t) < e^“^ g(t) . 

Thus 

*0^t<T ll®ll L^(fl) ^ hJ(«) 

* i llHllL2(a) PllL2(a)<“> > 

Now by the definition of U(x,0) 

and the desired result is obtained, that is 
0<t<T ll“-“llL2(n)'‘^> » I 

+ / lift (u-ujii^a,^, (t) dt ) 
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This implies that the error bounds are dependent on how well 
u and its first partial derivative with respect to time can 
be approximated by the arbitrary map u. For Hermite cubics, 
II u-u II l 2 (j^) is OCh"^) , and || u-u|| ^^1 is O(h^) . 

With further restrictions it can also be shown that the 
(fi) error is O(h^). 

In [3] J. E. Bendy, Jr. proves the following theorem 
which gives error estimates for the discrete time case. 

Theorem 2 . Let u be the solution of 

<1^ , v> + <a ( • , u) Vu, Vv> = <f(*,u),v> , veHQ(fi) 

u(x, 0) = Uq (x) X e 

u( • ,t) Hj(n) ^ , t e (0,T) 

and assume that ^ ^ 'J = max a (x, t , u) , (x, t)eHx (0 , t] . 

Let , where e M, and M satisfies 

inf ||u-u|| £ C h^”^ Hull [3]. Let W be the 

mapping from: (0,T] defined by 

<a(',u (u-W) ,V> + <u-W,V> = 0 for all VeM . 

2 

If ueC^(fi^(0,T]), and u , — 2" ^ (0,T;H {^) ) , r^2, 

3 1 

then for 1 < N < M and h and At sufficiently small 
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9 ||9e”“|| ^ At + q( ||e”|| ^ + ||e'' 1|| h 

+ 2A^At II (E^ - E^"^) II ^ ^ 

II u^-U^) I I = e^^^^(At)2 = c h^ , 

if 

llE^ll I + |lE°|| 2 + At |||^ (E^-Eq) II 2 < c(At)^ 

For the case of Hermite cubics, r = 4, and the conditions 
on M are satisfied. U is the approximate solution obtained 
in the Laplace modified alternating direction method. 

Using the program it can be shown that the order of the 
error is actually this by approximating the value of to where 
the error is O(h^). This approximation can be obtained in 
the following way; ' 

Let e^^ be the error for h = h^^, and B 2 be the error for 
h = h 2 / then 

e^ = C h“ e2 = C h2 

log(ej^/e 2 ) = to log(hj^/h 2 ) 

to = (log(ej^) - log(e 2 ) )/(log(h^) -log(h 2 )) 
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3 ; >9 



v;e obtain the 



] /2 

Using this method and choosing t = h 
results presented in Table 1. Table 2 shows similar results 
for the time error. 

To obtain the estimates in the tv;o tables the following 
problem was used for easily obtained analytic solution. 



Ut = - 2(e {(x-x^) + (y-y^) + ^(x-x^) (y-y^)) 

u(*,t) = 0 for all X e 



u(x,y,0) 



(x-x^) (y-y^) 



In this test 1 = -r v?as used in the discrete time 

4 

formulation . 



TABLE 1 



H 

0.2000 

0.1000 

THE VALUE 
THE VALUE 



DT 

0.0400 

0.0100 

OF W FOR 
OF VJ FOR 



L-2 NORM 

0.0003680 

0.0000286 

THE MAX NORM IS 
THE L-2 NORM IS 



MAX NORM 

0.0006815 

0.0000619 

.46006. 

.68510. 



TABLE 2 



H 

0.2000 

0.1000 

THE VALUE 
THE VALUE 



DT 

0.0400 

0.0100 



L-2 NORM 
0.0003680 
0.0000286 



I4AX NORl-1 
0.0006815 
0.0000619 



OF W FOR THE MAX NORM IS 1.73003. 
OF W FOR THE L-2 NORM IS 1.84255. 
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VI . TEST PROBLEMS 



The utility of the program is evident in several fields 
of engineering. The following test problems will illustrate 
its use in two of these. The problems are taken from the 
fields of fluid dynamics and heat transfer. 

A. RAYLEIGH'S PROBLEM FOR A CORNER 

In the first problem we consider the diffusion equation 
in two dimensions with the following boundary and initial 
conditions 

u^ = (1/Re) (u^^ + Uyy) , x,y > 0 
u ( X , y , 0 ) = 1 

u(0,y,t) = u(x,0,t) = 0 
u(x,y,t) ->-1 as X and y ~ 

where Re is the Reynold's number of the fluid under 
consideration. 

This is known as Rayleigh's problem for a corner, and 
the solution describes the impulsive motion of a right 
angled corner formed by two infinite flat plates and is 
used to infer the steady flow along the corner, with leading 
edge at t=0 . 



35 



For the purpose of numerical approximations, using a 
Reynold's number of 1000, x = 4 and y = 4 can be considered 
as infinity and the problem can be stated for use in the 
program as 

u. = (1/1000) (u + u ) 

t XX yy 

u(x,y,0) = 1 

u(0,y,t) = u(x,0,t) = 0 

u(4,y,t) = u(x,4,t) = 1 

and 

a(x,y,t,u) = 1/1000 
f (x,y, t,u,u^,u^) = 0 

A nonuniform grid with a fine mesh near the x = 0, y = 0 
corner is useful for an accurate description of the boundary 
layer behavior near the corner. Illustration (1) is a 
diagram of the grid used. 

The exact solution to this problem can be found analyt- 
ically to be 

u(x,y,t) = erf(X) erf (Y) 

where 

X = (x/2) (Re/t) 

Y = (y/2) (Re/t) 



36 



Illustration ( 1 ) 



The nature of the solution is that of a shock wave 



traveling through the rectangle. The error in the approxi- 
mation seems to be greatest on the wave and small even when 
the grid size is quite large. The error was also larger at 
early times near the boundary discontinuity than later in 
the process. Table 3 is a sample of the output of the 
program for this problem. The output points were taken near 
the discontinuity on the boundary since it is this area 
which is of greatest interest. 

For this problem X = .00025, At = .03 , and the 
largest Ax was .75. 

B. THE HEAT TRANSFER PROBLEM 

The equation for heat conductance for heterogeneous 
isotropic solids and frictionless incompressible fluids is 

pc |H = V (k u) + qQ- 

where the solution u(x,y,t) represents the temperature of 
the material. p is the density of the material, c is the 
heat capacitance, k is the thermoconductivity, and <gQ"' is 
an internal generation term. pc is a constant and can be 
combined with k and qQ'" to give 

= V- (k' u) + q"'/pc 

In this form k' is the thermal dif fusitivity . 
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For most materials the thermal dif fusitivity can be 
approximated by a linear function of the temperature. It 
is typically of the form 

k' = k + a(u-UQ) 

where k is the initial value of the thermal dif fusitivity 
and a is a small real number. 

The internal generation term is highly variable and 
dependent on the type of problem. A typical function 
encountered is of the form 



Analytic solutions for problems such as this are not 
available and it is in these cases that approximation 
techniques are invaluable. 

The problem for the program run will be as follows 



u^ = V ( (100 + .005 (u-Uq) ) Vu) + 

u(x,y,0) = sin(TTX) + sin(Tiy) 

u(0,y,t) = u(l,y,t) = ^ sin(Try) 

/_ 2 . 

u(x,0,t) = u(x,l,t) = e^ ' sin(iTx) 



a(x,y,t,u) = 100 t 
f (x,y , t,u,u^,u,^) = 



.005 (u- (sin(nx) +sin(7ry))) 
( u/ (sin ( TTx) +sin ( iry) ) ) 

G 
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with uniform grid, with Ax = Ay = h = 1/10 , time 
step At = .01 and X = 30 . 

Plots 1-8 show the general nature of the solution. 

Plot 1 shov?s the solution at an early time while it retains 
much of its initial form. The area near the boundary is 
beginning to become smaller. The internal generation term 
causes the oscillatory nature of the solution. After a 
longer period of time the solution begins to take on a 
steady state of zero. 
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PLOT (1) TIME = .05 
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/ 

/ 



PLOT (2) TIME .40 
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PLOT (3) TIME = .5 
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PLOT (4) TIME == .6 




PLOT (5) TIME = .7 
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/ 



PLOT (6) TIME = .9 
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PLOT (8) TIME = 2.55 




VII. CONCLUSIONS AND RECOMMENDATIONS 



Two major conclusions can be reached from the applica- 
tion of the program to the various problems. First the 
program is in general more efficient than the standard 
implementation of Galerkin methods which achieve the same 
degree of accuracy. It is also possible to achieve a higher 
degree of continuity and approximation in the solution using 
higher degree B-splines. It is possible that changes not 
evident to the author may be made to improve the efficiency 
of the program. Most of the computational effort in a given 
time step is in evaluating and VU^ at the gauss points; 
more efficient evaluations or placement of the quadrature 
points can possibly increase the speed of the program. It 
should be noted that the matrices [see (3. 8), (3. 9)] arising 
need be factored only once since they are time-independent. 
Second, optimal error bounds for the method have not been 
derived in the case that the Laplace modified forward 
difference method is used to obtain the second set of 
coefficients so that the centered difference method may be 
used, but this does not seem to have any effect on the 
solution in application. 

Using the program as a basis the Galerkin methods may 
be extended to a third dimension, using extension of the 
alternating direction to obtain it. The basic ideas used 
in the development of the program can be utilized to program 
the second order hyperbolic problem using the Galerkin methods. 
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APPENDIX A 



THE COMPUTER PROGRAM 

A. THE PROGRAM VARIABLES 

Alpham - the coefficients of the splines at time m. 
CX,CY,AX,AY - the matrices C„,C, .A^,A„ respectively 
described in Section III 
X,Y - the grid points on the x and y axes 
TKX,TKY - the knots system on the x and y axes, used 
for the de Boor routines 
GX/GY - the Gauss weights computed in WEIGHT 
EKCX,ETAY - the values of the basis functions at the 
quadrature points 

ALXR, ALXL , ALYB , ALYT - the coefficients for the projection 
of the boundary conditions at time m 
CPX,CPY - the quadrature points on the x and y axes 
NX, NY - the number of intervals on the x and y axes 
IAX,IAY - the dimension of the matrices defined above 
DT - At 
TIME - time m 

B. THE MAIN PROGRAM 

The main program calls various subroutines to accomplish 
the necessary procedures. Brief indications of the function 
of each subroutine will follow. The program begins by 
setting up the grid desired and its corresponding knots 
system and then computes the necessary matrices. It stores 
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parts of the matrices which will be changed but needed for 
future corrections. The first set of coefficients is 
obtained by projecting the initial conditions into the 
space. The second set is computed using a forward differ- 
ence Laplace modified step. The program then steps forward 
in time using the central difference Laplace modified 
method . 

C. MATCOM 

MATCOM computes the matrices for the program and stores 
them in band symmetric storage used by IMSL library routines. 
It also stores the quadrature points and the values of the 
basis factions at the quadrature points. 

D. FMATCO 

FMATCO projects the initial conditions into the space. 

E. \VEIGHT 

WEIGHT computes and stores the product of the quadrature 
weights and Ax and Ay for the appropriate interval. 

F. GRID 

GRID computes the grid points. 

G. KNOTS 

KNOTS computes the knot system for the de Boor routines. 

H . UEVAL 

UEVAL computes the value of the function at specific 
points (because of inefficiency its use should be limited) . 
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I. MULTI P 

MULTIP computes (CX CY) in manner described in 

Section III. 

J. SOLVER 

SOLVER solves system of the form (CX I) (I CY) e^ = 3 
in the manner described in Section III. 

K. ANEVAL 

ANEVAL evaluates the function and its first partial 
derivatives at the quadrature points. 

L. RTHDSD 

RTHDSD computes the right hand side of equation (3.8). 

M. ALXYBC 

ALXYBC projects the boundary conditions into the space. 

N. FUNCTION F 

Function F defines the necessary functions for the 
problem being solved. 

O. THE DE BOOR SUBROUTINES 

The de Boor subroutines perform various operations to 
obtain B-splines [2]. 
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APPENDIX B 



THE I? PROJECTION OF FUNCTIONS 
2 . . 

The L projection of the initial conditions and other 
functions involved in the program were accomplished using 
Gauss quadrature. Appendix A describes several subroutines 
in which this is done. To achieve the desired accuracy it 
was necessary to use nine quadrature points for each interval 
In the subroutine \'JEIGHT, the Gauss weights are combined with 
the value of Ax and Ay for their respective intervals and 
saved. This is necessary in the case of non-uniform grids. 
The quadrature points are stored for future use in the 
subroutine MATCOM. 

To project a function into the space, we must take its 
inner produce with each of the basis functions. This means 
we must know the values of the basis f\inctions at the 
quadrature points. These are computed using the de Boor 
routines and saved in MATCOM. In the case of one dimensional 
B-splines there are only four non-zero basis functions on 
each interval and each basis function is non-zero on no more 
than two intervals. Thus on each interval the quadrature is 
taken with the functions and each of the four non-zero basis 
functions . 

+" Vj Vi 

For example, on the I x-interval and the J y-interval 
we evaluate the function at the nine quadrature points and 
take its product with the value of the basis functions at 
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the point and the appropriate Gauss weights, and then sum 
over the nine points in the interval for each of the four 
basis functions. This procedure is continued by moving to 
the next interval, and using the appropriate basis functions 
for that interval. 
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